function CosGibbsMH(simp,Qd,m,s,Bin_DepVar,Score_assessor,Sc_Undes,iternum,model)
    %implements the MCMC simulation using the Metropolis-Hastings (M-H) algorithm to generate a sequence of random samples from the posterior distribution.  
    %"Cos-": addresses the stratified sampling problem using the approach of Cosslett (1993)
    %"-Gibbs-": using Gibbs sampling of parameters (mu, sigma) for covariate X  [See Appendix B in Lieli and Springborn]
    %"-MH": using the Metropolis-Hastings algorithm

%Inputs to CosGibbsMH()
    %Qd: Assumed base rate of weeds (unconditional probability that Y=1).  Note: referred to in manuscript as "tau".
    %simp: Model for the density of the covariate sample, f(x;lambda) in Equation (13) of Lieli and Springborn.  [Also see Appendix B in Lieli and Springborn]
        %Reported results in Lieli and Springborn are for simp=1 (mean and variance of x set to their maximum likelihood estimates)
        %See "CosGibbsMH.m" for detail on all four variations of increasing complexity.
    %iternum: number for a particular random subset of the data.  Only used if running model on random subsets of the data.
    %m: chain index (will simulate several different Markov chains).  
    %s: iterations in each chain
    %Bin_DepVar,Score_assessor,Sc_Undes: data
    %model: link function (e.g. logit, cauchit).
    

%Results:
    %Output is saved to a particular folder structure which should be created ahead of time:
        %results\[data used]\Q*\simp*\[model]\[filename].mat
            %[data used] may be either "mcmc_univ_subsample" if using subsamples of data, or "mcmc_univariate" if using the whole data set
            %full titles for "Q*" and "simp*" will be given by the particular base rate (Q) and covariate distribution simplification (simp) being run
            %[model] is given by the link function.
            %[filename] is given by several variables (see below) including the chain.             
    
%************************************

 
%% Set simplification (1:4) corresponding to complexity of model for f(x|y,muy,sigy]).
    %1:  [Nonrandom] Set mu's & sig's to modal levels
    %2:  [mu's random; sig's fixed and shared]  f(x|y=0) and f(x|y=1) have same sig's (fixed at mode) but different mu's random
    %3:  [mu's random; sig's random and shared] Variance of x|y is random but the same for y=1 and y=0, i.e. draw of sig is shared.
    %4:  [mu's random; sig's random; not shared] This is the fully random model.
simplification=simp;%

%% Likelihood functions
%Depending on the link function (logit or cauchit) specify function for calculating likehlihood component of posterior.
if strcmp(model,'logit')
    logitp = @(b0,b1,x) 1./(1+exp(-(b0+b1.*x)));%Need this form since typical form returns NaN for high b1 (e.g. b1=1000): exp(b0+b1.*x)./(1+exp(b0+b1.*x));  %
    fygvnXB = @(b0,b1,x,y) (logitp(b0,b1,x).^y).*((1-logitp(b0,b1,x)).^(1-y)); %f(y|x,B), logit link
elseif strcmp(model,'cauchit')
    cauchitp= @(b0,b1,x) .5 + atan(b0+b1.*x)/pi; %ezplot(@(x) cauchitp(.14,.3,x),[-10,10])
    fygvnXB = @(b0,b1,x,y) (cauchitp(b0,b1,x).^y).*((1-cauchitp(b0,b1,x)).^(1-y)); 
end

%standard "raw" likelihood in absence of stratified samp. issue.
rawlike = @(b0,b1,x,y) exp(sum(log(fygvnXB(b0,b1,x,y))));  


%% Proposal density for beta1: 
propdraw=@(mubeta1,X,weta,k) mubeta1 + weta*normrnd(0,1,k,1);

%% Function for imposing the constraint on Qd (tau).  [See Lieli and Springborn Equation (14)]
%Zero function for satisfying the constraint (square the difference to facilitate zero-finding).  
%Note: fygvnXB(b0,b1,x,y=1)=logitp(b0,b1,x)
lb=-10^3; ub=10^3;  %upper and lower bounds for the integration.  May be unstable outside these bounds.  Bounds over x provide sufficient precision.
if strcmp(model,'logit')
    constr=@(b0,b1,mu0,sig0,mu1,sig1,lb,ub,Qd) (Qd - (quadgk(@(x) logitp(b0,b1,x).*normpdf(x,mu0,sig0), lb,ub))./(1 - quadgk(@(x) logitp(b0,b1,x).*(normpdf(x,mu1,sig1)-normpdf(x,mu0,sig0)), lb,ub))).^2; %
elseif strcmp(model,'cauchit')
    constr=@(b0,b1,mu0,sig0,mu1,sig1,lb,ub,Qd) (Qd - (quadgk(@(x) cauchitp(b0,b1,x).*normpdf(x,mu0,sig0), lb,ub))./(1 - quadgk(@(x) cauchitp(b0,b1,x).*(normpdf(x,mu1,sig1)-normpdf(x,mu0,sig0)), lb,ub))).^2; %
end    
    
 
%% The prior over theta.
%Using a "data augmentation prior" (DAP).   
%Since using a uniform prior, these vectors will be empty. 
X_dap=[];%
if strcmp(model,'logit')
    y_dap=[];%
elseif strcmp(model,'cauchit')
    y_dap=[]; %
end
w_dap=[]; %weights on prior. 


    
%% MCMC analysis
%Set constants and arrays
k=1; %number of parameters in M-H step
betatarget=0.4;   %Target M-H step acceptance rate: Roberts, Gelman and Gilks show approx .45 is good for one-dimensional problems.
wetatol   =1.1;   %Tuning parameter for dispersion of proposal density
weta=.75;         %starting point for dispersion parameter in M-H step. 
                  
%% setup for using just a sample of observations to enable out of sample calculations...
if isempty(iternum)
    subsamp=0;
else
    %subsampmat is a 370x1000 matrix (see runBayesEst.m) where each vector
    %contains the numbers 1:370 in random order.  
    subsamp=1;
    temp=load('D:\Bayesian_Analysis\subsampmat.mat');  %generated by runBayesEst.m
    ssvec=temp.ssmat(1:250,iternum);  %Get record numbers for 250 randomly chosen observations.
end
    
%%Specify data to use: either whole or sub-sample of the data
if subsamp==0 %Use all the data 
    y=Bin_DepVar; 
    X=[Score_assessor]; %no need to include ones given form of logitp
elseif subsamp==1  %Use subsample of the data
    y=Bin_DepVar(ssvec); 
    X=[Score_assessor(ssvec)]; %no need to include ones given form of logitp
end

%% Set params according to the data:
N=length(y); %Number of observations in training data.
n1=sum(y); n0=sum(1-y);
ybar=mean(y); 
H=ybar;  %Sampling intensity of observations for which y=1

%Implement "data augmenting" form of prior. 
%Here the prior is assumed uniform so _dap vectors are blank.
X_aug=[X;X_dap];
y_aug=[y; y_dap.*w_dap];

%%Set empty matricies to store results
beta1s=zeros(s,k); %array for the accepted simulation draws
beta0s=zeros(s,k); %array for constant term
% stats will have prefix v. 
vb1mean =zeros(s,k);  
vb1std  =zeros(s,k);
vb0mean =zeros(s,k); % stats will have prefix v.  
vb0std  =zeros(s,k);
petas =zeros(s,k);
b0ofpetas=zeros(s,k);
wetas=zeros(s,k);
hbeta=0; %counter for # of accepted jumps

%% Draw Gibbs samples from distribution parameters for covariate, X.
%Note: suffix "0" denotes terms that pretain to f(x|y=0), similarly for suffix "1".
%Structure:  Sigma0^2 distributed scale-Inv-Chi2(v0,s0^2), mu0 distributed N(mean(X(y==0)), sigma0^2).  Similarly for y=1.
%Random draw var from Scale-Inv-Chi2(v,s^2): draw Z from Chi2(v), then var=v*s^2/Z, and sigrnd0=var.^(.5)
v0=n0-1; s20=std(X(y==0))^2; 
v1=n1-1; s21=std(X(y==1))^2;

%Set mus and sigmas based on approach specified in "simplification":
    %1:  [Nonrandom] Set mu & sig to modal levels
    %2:  [mu's random; sig's fixed and shared]  f(x|y=0) and f(x|y=1) have same sig's (fixed at mode) but different mu's random
    %3:  [mu's random; sig's random and shared] Variance of x|y is random but the same for y=1 and y=0, i.e. draw of sig is shared.
    %4:  [mu's random; sig's random; not shared]
    
if simplification==1 % set f(x) parameters to modal levels: 
    sig0rnd=sqrt(((v0*s20)/(v0+2)))*ones(s,1); mu0rnd=mean(X(y==0))*ones(s,1);
    sig1rnd=sqrt(((v1*s21)/(v1+2)))*ones(s,1); mu1rnd=mean(X(y==1))*ones(s,1);
elseif simplification==2 || simplification==3 %simplified draws on x.  Assumes dist of x for y=1 and y=0 is the same accept for a shift in the mean.  
    mudraw=randn(s,1);
    Xdemean((y==0))=X(y==0) - mean(X(y==0));Xdemean((y==1))=X(y==1) - mean(X(y==1));
    v = N-1; s2=std(Xdemean)^2;
    if simplification==2 %Variance is fixed/known
        sig0rnd=sqrt(((v*s2)/(v+2)))*ones(s,1); 
        sig1rnd=sig0rnd;
    elseif simplification==3 
        Z = chi2rnd(v,s,1); sigrnd=sqrt((v*s2)./Z); sig0rnd=sigrnd; sig1rnd=sigrnd;
    end
    mu0rnd=mean(X(y==0))+ sig0rnd.*mudraw;
    mu1rnd=mean(X(y==1))+ sig1rnd.*mudraw;
elseif simplification==4
    Z0=chi2rnd(v0,s,1); sig0rnd=sqrt((v0*s20)./Z0); mu0rnd=mean(X(y==0))+ sig0rnd.*randn(s,1);   %
    Z1=chi2rnd(v1,s,1); sig1rnd=sqrt((v1*s21)./Z1); mu1rnd=mean(X(y==1))+ sig1rnd.*randn(s,1);   %
end

%Need to exclude cases in which draw of mu0rnd > mu1rnd (which is infeasible).   
ind=find(mu0rnd>mu1rnd); length(ind)
while ~isempty(ind)
    mu0rnd(ind)=mean(X(y==0))+ sig0rnd(ind).*randn(length(ind),1);
    mu1rnd(ind)=mean(X(y==1))+ sig1rnd(ind).*randn(length(ind),1);
    ind=find(mu0rnd>mu1rnd);
end

 
%% Initialize chain: Set starting values for the 11 chains
if strcmp(model,'logit')
    %Starting values for beta parameters: Using maximum likelihood estimates (from runmle.m).  
    %Begin with MLE estimates uncorrected for stratified sample:
    bmleuncorr=[.1457 .29450]; 
    %For logit, correction for stratified sample can be implemented after the fact 
    %Using the King and Zeng correction on constant takes care of shifts in Q:
    bmle=[(bmleuncorr(1) - log(((1-Qd)/Qd)*(ybar/(1-ybar)))) bmleuncorr(2)]; %King and Zeng prior correction: B0-ln(stuff)  pg 1415, stat in medicine 2002.
    %NOTE: the King and Zeng correction is necessary for the constant term (b0) but not the slope term (b1).
elseif strcmp(model,'cauchit') %For Cauchit, King & Zeng correction not approp. so use Manski-Lerman weights in runmle.m to get baseline est.
    if Qd==0.02
        bmle=[-31.2032  2.4559]; 
    elseif Qd==0.05
        bmle=[-12.2617  0.9860]; 
    end
end

%calculate the starting point for b1 for M=11 chains at the specified quantiles.
b1start = norminv([.01 .05 .1 .2 .35 .5 .65 .8 .9 .95 .99],bmle(2),.33); %.33 was the sd observed for Q=0.02 starting at mean.  
%initialize chain:
beta1s(1)=b1start(m); %for the m-th chain
beta1=beta1s(1); 

%Find beta0 that satisfies the constraint on Q (tau) given the corrected b0 from above (bmle(1)) as a starting point.
[beta0s(1),fval(1),exitflag(1)] = fminsearch(@(b0) constr(b0,beta1,mu0rnd(1),sig0rnd(1),mu1rnd(1),sig1rnd(1),lb,ub,Qd),bmle(1)); beta0=beta0s(1);
wetas(1)=weta;
arbeta=0;
    
%% Run the MCMC 
for t=2:s
    %%Gibbs sampling step: mu and sigma as given by: mu0rnd(t),sig0rnd(t),mu1rnd(t),sig1rnd(t)... previously drawn as independent.
    %Metropolis-Hastings step:
    peta=propdraw(beta1s(t-1),X,weta,k);  %draw a proposal for beta1 from the proposal density
    %Given Qd and peta, solve for constant term b0:
    [b0ofpeta,fval(t),exitflag(t)] = fminsearch(@(b0) constr(b0,peta,mu0rnd(t),sig0rnd(t),mu1rnd(t),sig1rnd(t),lb,ub,Qd),beta0s(t-1));

    %For the acceptance probability, I need the current beta0 & beta1 in the chain.  But since the parameters of f(x)
    %have changed and Qd is fixed, I need to recalculate beta0(beta1, mu(t), sig(t)).
    %If I treat mu & sig as fixed this step will not be needed.
    if simp==1 %If simp==1, i.e. if mu & sigma fixed then b0 corresponding to beta1 doesn't change.
        b0ofbeta1=beta0s(t-1);
    elseif simp~=1
        [b0ofbeta1,fval2(t),exitflag2(t)] = fminsearch(@(b0) constr(b0,beta1,mu0rnd(t),sig0rnd(t),mu1rnd(t),sig1rnd(t),lb,ub,Qd),beta0s(t-1)); 
    end
    
    %%Accept or reject the proposed jump:        
    %Calculate log version of acceptance probability:        
    logratio= sum(log(fygvnXB(b0ofpeta   ,peta       ,X_aug,y_aug)))-sum(log(fygvnXB(b0ofbeta1,beta1,X_aug,y_aug)));% ...  %Note: portions of the posterior not a function of theta were excluded since they will cancel in the ratio
            %- log(propdens(peta       ,beta1s(t-1),X,weta)) + log(propdens(beta1s(t-1),peta       ,X,weta));  %Note: proposal density portion of acceptance ratio not needed since it's symetric and cancels.
            %above collapses to Metropolis algorithm...
    if binornd(1,exp(min(logratio,0)))==1  %Jump is accepted
        beta1=peta;     %update current state of chain if jump accepted
        beta0=b0ofpeta; %update current state of chain if jump accepted
        hbeta=hbeta+1;  %a count of accepted jumps
    end
    petas(t,:)    =peta';
    b0ofpetas(t,:)=b0ofpeta';

    beta1s(t,:)=beta1';
    beta0s(t,:)=beta0';

    wetas(t,:)  =weta';
    arbetas(t,:)=arbeta';

    arbeta=hbeta/t;     %rate of jump acceptance 

    lbw=0; ubw=100; %Q=0.02, fx variable.
    intw=10; %every intw rounds, revisit weta
    if t/intw==round(t/intw)
        dev=(arbeta-betatarget)/betatarget;
        weta=weta*exp(dev);
    end
    %calculate stats
    vb1mean(t)= mean(beta1s(1:t)) ; % stats will have prefix v.  
    vb1std  =zeros(s,k);
    vb0mean(t) =mean(beta0s(1:t)); % stats will have prefix v.  
    vb0std  =zeros(s,k);

    %Monitoring output.  Uncomment code to see figure.  
    inter=100;
    if t/inter == round(t/inter) %every inter-th iteration update plot
         display(['Iter#=' num2str(iternum) ', Q=' num2str(Qd) ', simp=' num2str(simp) ', m=' num2str(m) ', t=' num2str(t) '/' num2str(s)])
%             figure(100); clf; %use clear figure (clf) to avoid overlapping plots and killing memory.
%             subplot(4,1,1); ylabel('b0'); stairs(2:t,vb0mean(2:t),'k'); hold on; stairs(1:t,beta0s(1:t),'b'); 
%             subplot(4,1,2); ylabel('b1'); stairs(2:t,vb1mean(2:t),'k'); hold on; stairs(1:t,beta1s(1:t),'b'); 
%             subplot(4,1,3); ylabel('Acceptance Rate'); stairs(2:t,arbetas(2:t),'k');
%             subplot(4,1,4); ylabel('weta'); plot(2:t,wetas(2:t),'k-');
%             if t==inter | t==inter*2; pause(.1); end
    end
end  %end loop on t
    
    
if subsamp==0
    folder='mcmc_univariate';
    itertext=[];
elseif subsamp==1
    folder='mcmc_univ_subsample'
    itertext=['_iter' num2str(iternum)];
end
dest=['D:\Bayesian_Analysis\results\' folder '\Q' num2str(Qd) '\simp' num2str(simplification) '\' model '\Q' num2str(Qd) '_simp' num2str(simplification) '_m' num2str(m) '_s' num2str(s) itertext '.mat']
save(dest);
  
    
return



    